Cluster Monte Carlo and numerical mean field analysis for the water liquid— liquid 

phase transition 



OS 

o 
o 

(N 

G 

c3 



Marco G. Mazza, 1 Kevin Stokely, 1 Elena Strekalova, 1 H. Eugene Stanley, 1 and Giancarlo Franzese 2 

1 Center for Polymer Studies and Department of Physics, 
Boston University, Boston, Massachusetts 02215, USA 
2 Departament de Fisica Fonamental, Universitat de Barcelona, Diagonal 647, 08028 Barcelona, Spain 

By the Wolff's cluster Monte Carlo simulations and numerical minimization within a mean field 
approach, we study the low temperature phase diagram of water, adopting a cell model that repro- 
duces the known properties of water in its fluid phases. Both methods allows us to study the water 
thermodynamic behavior at temperatures where other numerical approaches -both Monte Carlo 
and molecular dynamics- are seriously hampered by the large increase of the correlation times. The 
cluster algorithm also allows us to emphasize that the liquid-liquid phase transition corresponds to 
the percolation transition of tetrahedrally ordered water molecules. 

PACS numbers: 61.20.Ja, 61.20.Gy 
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INTRODUCTION 

Water is possibly the most important liquid for life 
1] and, at the same time, is a very peculiar liquid Q. 
In the stable liquid regime its thermodynamic response 
functions behave qualitatively differently than a typi- 
cal liquid. The isothermal compressibility Kj-, for ex- 
ample, has a minimum as a function of temperature at 
T = 46 °C, while for a typical liquid monotonically 
decreases upon cooling. Water's anomalies become even 
more pronounced as the system is cooled below the melt- 
ing point and enters the metastable supercooled regime 
1- 

Different hypothesis have been proposed to rational- 
ize the anomalies of water All these interpretations, 
but one, predict the existence of a liquid-liquid phase 
transition in the supercooled state, consistent with the 
experiments to date [4( and supported by different mod- 
els 0. 

To discriminate among the different interpretations, 
many experiments have been performed [5[. However, 
the freezing in the temperature-range of interest can be 
avoided only for water in confined geometries or on the 
surface of macromolecules 0, @. Since experiments in 
the supercooled region are difficult to perform, numeri- 
cal simulations have played an important role in recent 
years to help interpret the data. However, also the sim- 
ulations at very low temperature T are hampered by the 
glassy dynamics of the empirical models of water 0, Q . 
For these reasons is important to implement more effi- 
cient numerical simulations for simple models, able to 
capture the fundamental physics of water but also less 
computationally expensive. Here we introduce the im- 
plementation of a Wolff's cluster algorithm 9] for the 
Monte Carlo (MC) simulations of a cell model for wa- 
ter The model is able to reproduce all the differ- 
ent scenarios proposed to interpret the behavior of wa- 
ter llj and has been analyzed (i) with mean field (MF) 



[HHHm, (ii) with Metropolis MC simulations @, EH 
and (iii) with Wang-Landau MC density of state algo- 
rithm [l5j |. Recent Metropolis MC simulations [|[ have 
shown that very large times are needed to equilibrate 
the system as T — > 0, as a consequence of the onset of 
the glassy dynamics. The implementation of the Wolff's 
clusters MC dynamics, presented here, allows us to (i) 
drastically reduce the equilibration times of the model at 
very low T and (ii) give a geometrical characterization of 
the regions of correlated water molecules (clusters) at low 
T and show that the liquid-liquid phase transition can 
be interpreted as a percolation transition of the tetrahe- 
drally ordered clusters. 

THE MODEL 

The system consists of N particles distributed within a 
volume V in d dimensions. The volume is divided into iV 
cells of volume Vi with i £ [1, N], For sake of simplicity, 
these cells are chosen of the same size, Vi ~ V/N, but the 
generalization to the case in which the volume can change 
without changes in the topology of the nearest-neighbor 
(n.n.) is straightforward. By definition, vi > vq, where 
vq is the molecule hard-core volume. Each cell has a 
variable = for a gas-like or rii = 1 for a liquid-like 
cell. We partition the total volume in a way such that 
each cell has at least four n.n. cells, e.g. as in a cubic 
lattice in 3d or a square lattice in 2d. Periodic boundary 
conditions are used to limit finite-size effects. 

The system is described by the Hamiltonian [1(3] 

-Ja^rii (1) 

i (k,l)i 

where e > is the strength of the van der Waals attrac- 
tion, J > accounts for the hydrogen bond energy, with 



2 




FIG. 1: A pictorial representation of five water molecules in 
3d. Two hydrogen bonds (grey links) connect the hydrogens 
(in blue) of the central molecule with the lone electrons (small 
gray lines) of two nearest neighbor (n.n.) molecules. A bond 
index (arm) with q = 6 possible values is associated to each 
hydrogen and lone electron, giving rise to q 4 possible orien- 
tational states for each molecule. A hydrogen bond can be 
formed only if the two facing arms of the n.n. molecules are 
in the same state. Arms on the same molecule interact among 
themselves to mimic the O-O-O interaction that drives the 
molecules toward a tetrahedral local structure. 

four (Potts) variables cry = 1, . . . , q representing bond in- 
dices of molecule i with respect to the four n.n. molecules 
j, 3a,b — 1 if a — b and 5 a ,b = otherwise, and {i,j) de- 
notes that i and j are n.n. The model does not assume a 
privileged state for bond formation. Any time two facing 
bond indices (arms) are in the same (Potts) state, a bond 
is formed. The third term represents an intramolecular 
(IM ) interaction accounting for the 0-0-0 correlation 
[161 ] . locally driving the molecules toward a tetrahedral 
configuration. When the bond indices of a molecule are 
in the same state, the energy is decreased by an amount 
J„ ^ and we associate this local ordered configuration 
to a local tetrahedral arrangement [17] ■ The notation 
(k, l)i indicates one of the six different pairs of the four 
bond indices of molecule i (Figfl]). 

Experiments show that the formation of a hydrogen 
bond leads to a local volume expansion Q. Thus in our 
system the total volume is 

V = Nv + N HB v HB , (2) 

where 

N H b = ^ n i n J^tj^n ( 3 ) 
<i,j> 



constant specific volume increase due to the hydrogen 
bond formation. 

MEAN-FIELD ANALYSIS 

In the mean-field (MF) analysis the macrostate of 
the system in equilibrium at constant pressure P and 
temperature T (NPT ensemble) may be determined by 
a minimization of the Gibbs free energy per molecule, 
g = {(Ji?) + PV - TS)/N W , where 

N w = J2 n * ( 4 ) 

i 

is the total number of liquid-like cells, and S = S n + S a 
is the sum of the entropy S n over the variables rij and 
the entropy S a over the variables cry . 

A MF approach consists of writing g explicitly using 
the approximations 

n i n o — ► 2Nn2 (5) 

<ij> 

^ niUjS^^i — ► 2Nn 2 p a (6) 

<ij> 

"W.^i — * 6Nnp a (7) 

i (k,l)i 

where n — N w /N is the average of n,, and p a is the 
probability that two adjacent bond indices cry are in the 
appropriate state to form a hydrogen bond. 
Therefore, in this approximation we can write 

V = Nv + 2Nn 2 PtT v H B, (8) 
(JF) = -2 [en +{Jn + 3 J a ) <p a ] nN. (9) 

The probability p,y, properly defined as the thermody- 
namic average over the whole system, is approximated 
as the average over two neighboring molecules, under the 
effect of the mean-field h of the surrounding molecules 

^ = <<WA- ( 10 ) 

The ground state of the system consists of all N vari- 
ables n,- = 1, and all cry in the same state. At low 
temperatures, the symmetry will remain broken, with 
the majority of the cry in the preferred state. We as- 
sociate this preferred state to the tetrahedral order of 
the molecules and define m a as the density of the bond 
indices in the tetrahedral state, with < m a < 1. There- 
fore, the number density n„ of bond indices cry is in the 
tetrahedral state is 

1 + (g - l)m a 
n a = . (11) 

Since an appropriate form for h is (loj 



is the total number of hydrogen bonds, and vrb is the 



(12) 
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we obtain that ^ <h< 3J a . 

The MF expressions for the entropies S n of the N vari- 
ables rii, and S a of the ANn variables then [H 

ff n = W(n log(n) + (1 - n) log(l - n)) (13) 



S'cr = -fcs4-/Vn[n CT log(n CT )+ 

(l-n CT )log(l-n CT )+log(g-l)], 

where ks is the Boltzmann constant. 
Equating 



(14) 



Po 



(1 - n a f 



(15) 



with the approximate expression in Eq. (jlOp . allows for 
solution of n a , and hence g, in terms of the order param- 
eter m a and n. 

By minimizing numerically the MF expression of g 
with respect to n and m a , wc find the equilibrium values 
n( e<? ) and and, with Eqs. (@|) and ((2]), we calcu- 

late the density p at any (T, P) and the full equation 
of state. An example of minimization of g is presented 
in Fig. [5] where, for the model's parameters J/e = 0.5, 
J CT /e = 0.05, vhb/vq — 0.5, q — 6, a discontinuity 
in m^ 69 -* is observed for Pvg/e > 0.8. As discussed in 
Ref.s | this discontinuity corresponds to a first or- 

der phase transition between two liquid phases with dif- 
ferent degree of tetrahedral order and, as a consequence, 
different density. The higher P at which the change 
in is continuous, corresponds to the pressure of a 

liquid-liquid critical point (LLCP). The occurrence of the 
LLCP is consistent with one of the possible interpreta- 
tions of the anomalies of water, as discussed in Ref. [12] • 
However, for different choices of parameters, the model 
reproduces also the other proposed scenarios 11 1 . 



THE SIMULATION WITH THE WOLFF'S 
CLUSTERS MONTE CARLO ALGORITHM 

To perform MC simulations in the NPT ensemble, we 
consider a modified version of the model in which we 
allow for continuous volume fluctuations. To this goal, 
(i) we assume that the system is homogeneous with all 
the variables rii set to 1 and all the cells with volume 
v = V/N; (ii) we consider that V = Vmc + NhbVhb, 
where Vmc Nvq is a dynamical variable allowed to 
fluctuate in the simulations; (iii) we replace the first (van 
der Waals) term of the Hamiltonian in Eq. fl} with a 
Lennard-Jones potential with attractive energy e > J 
and truncated at the hard-core distance 



U w (r) 



if r < r , 
if r > r . 



(16) 
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FIG. 2: Numerical minimization of the molar Gibbs free en- 
ergy g in the mean field approach. The model's parameters 
are J/e = 0.5, Jo-/e — 0.05, vhb/vo = 0.5 and q = 6. In each 
panel we present g (dashed lines) calculated at constant P 
and different values of T. The thick line crossing the dashed 
lines connects the minima m^ 9 ' of g at different T. Upper 
panel: Pv /e = 0.7, for T going from k B T/e = 0.06 (top) to 
k B T/e = 0.08 (bottom). Middle panel: Pv /e = 0.8, for T 
going from k B T/e = 0.05 (top) to k B T/e = 0.07 (bottom). 
Lower panel: Pvo/e — 0.9, for T going from k B T/e = 0.04 
(top) to k B T/e = 0.06 (bottom). In each panel dashed lines 

are separated by k B ST/e — 0.001. In all the panels m^ 9 ' 
increases when T decreases, being (marking the absence of 
tetrahedral order) at the higher temperatures and ~ 0.9 (high 
tetrahedral order) at the lowest temperature. By changing T, 



changes in a continuous way for Pvo/e — 0.7 and O.i 



but discontinuous for Pvo/e — 0.9 and higher P. 
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where tq = (vo) 1 ^; the distance between two n.n. 
molecules is (V/N) 1 ' d , and the distance r between two 
generic molecules is the Cartesian distance between the 
center of the cells in which they are included. 

The simplification (i) could be removed, by allowing 
the cells to assume different volumes Vi and keeping fixed 
the number of possible n.n. cells. However, the results 
of the model under the simplification (i) compares well 
with experiments [12]. Furthermore, the simplification 
(i) allows to drastically reduce the computational cost 
of the evaluation of the Uw{r) term from N(N — 1) to 
N—l operations. The changes (i)-(iii) modify the model 
used for the mean field analysis and allow off-lattice MC 
simulations for a cell model in which the topology of the 
molecules (i.e. the number of n.n.) is preserved. The 
comparison of the mean field results with the MC simu- 
lations show that these changes do not modify the physics 
of the system. 

We perform MC simulations with N = 2500 and 
N = 10000 molecules, each with four n.n. molecules, 
at constant P and T, in 2d, and with the same parame- 
ters used for the mean field analysis. To each molecules 
we associate a cell on a square lattice. The Wolff's algo- 
rithm is based on the definition of a cluster of variables 
chosen in such a way to be thermodynamically correlated 



191 ] . To define the Wolff's cluster, a bond index (arm) 
of a molecule is randomly selected; this is the initial el- 
ement of a stack. The cluster is grown by first checking 
the remaining arms of the same initial molecule: if they 
are in the same Potts state, then they are added to the 
stack with probability p sa me = min [1, 1 — exp(— /3J a )] 
[§], where (3 = (fc^T) -1 . This choice for the proba- 
bility Psame depends on the interaction ,J a between two 
arms on the same molecule and guarantees that the 
connected arms are thermodynamically correlated [l9j |. 
Next, the arm of a new molecule, facing the initially 
chosen arm, is considered. To guarantee that connected 
facing arms correspond to thermodynamically correlated 
variables, is necessary to link them with the probabil- 
ity Pfacing = min [1,1- cxp(-/3J')] where J' = J-Pv HB 
is the P-dependent effective coupling between two facing 
arms as results from the enthalpy J4? + PV of the sys- 
tem. It is important to note that J' can be positive or 
negative depending on P. If J' > and the two facing 
arms are in the same state, then the new arm is added 
to the stack with probability pf ac ing; if J' < and the 
two facing arms are in different states, then the new arm 
is added with probability j»f ; 



acing 



Only after every 
possible direction of growth for the cluster has been con- 
sidered the values of the arms are changed in a stochastic 
way; again we need to consider two cases: (i) if J' > 0, 
all arms are set to the same new value 



(cr old + </>) mod q 



(17) 



the same random constant <fi £ [1, . . . q] 

C = (a? d + </>) mod q. 



(18) 



In order to implement a constant P ensemble we let 
the volume fluctuate. A small increment Ar/r^ = 0.01 
is chosen with uniform random probability and added 
to the current radius of a cell. The change in volume 
AV = V ncw - V old and van der Waals energy AE W 
is computed and the move is accepted with probability 
min (1, exp [—(3 (AE W + PAV - TAS 1 )]), where AS = 
-Nk B ln(V new /V° ld ) is the entropic contribution. 



MONTE CARLO CORRELATION TIMES 

The cluster MC algorithm described in the previous 
section turns out to be very efficient at low T, allow- 
ing to study the thermodynamics of deeply supercooled 
water with quite intriguing results [2l| . To estimate the 
efficiency of the cluster MC dynamics with respect to the 
standard Metropolis MC dynamics, we evaluate in both 
dynamics, and compare, the autocorrelation function of 
the average magnetization per site Mi = \ Ylj a ij > where 
the sum is over the four bonding arms of molecule i. 



(M a fa+t)M t fa)) - (Mtf 
(Mf) - (Af,) 2 



(19) 



where <fi is a random integer between 1 and q; (ii) if J' < 
0, the state of every single arm is changed (rotated) by 



For sake of simplicity, we define the MC dynamics au- 
tocorrelation time r as the time, measured in MC steps, 
when Cm(t) = 1/e. Here we define a MC step as AN up- 
dates of the bond indices followed by a volume update, 
i.e. as 47V + 1 steps of the algorithm. 

In Fig. [3] we show a comparison of Cm(^) for the 
Metropolis and Wolff algorithm implementations of this 
model for a system with N — 50 x 50, at three tempera- 
tures along an isobar below the LLCP, and approaching 
the line of the maximum, but finite, correlation length, 
also known as Widom line Tw{P) [12fl . In the top panel, 
at T > T W {P) (k B T/e = 0.11, Pv /e = 0.6), we find 
a correlation time for the Wolff's cluster MC dynamics 
r\v ~ 3 x 10 3 , and for the Metropolis dynamics tm ~ 10 6 . 
In the middle panel, at T > T W (P) {k B T/e = 0.09, 
Pvq/c = 0.6) the difference between the two correlation 
times is larger: r w w 2.5 x 10 3 , t m ~ 3 x 10 6 . The bot- 
tom panel, at T ~ T W (P) {k B T/e = 0.06, Pv /e = 0.6) 
shows tw ~ 3.7 x 10 2 , while tm is beyond the accessible 
time window (tm > 10 7 ). 

Since as T — » the system enters a glassy state [1] , the 
efficiency tm/tw grows at lower T allowing the evalua- 
tion of thermodynamics averages even at T <C Tc 21 1. 
In particular, the cluster MC algorithm turns out to be 
very efficient when approaching the Widom line in the 
vicinity of the LLCP, with an efficiency of the order of 
10 4 . We plan to analyze in a systematic way how the 
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FIG. 3: Comparison of the autocorrelation function Cm(£) for 
the Metropolis (circles) and Wolff (squares) implementation of 
the present model. We show the temperatures fel/t = 0.1 1 
(top panel), k B T/e = 0.09 (middle panel), k B T/e = 0.06 
(bottom panel), along the isobar Pvo/e = 0.6 close to the 
LLCP for N = 50 x 50. 



efficiency tm/t\v grows on approaching the LLCP. This 
result is well known for the standard liquid-gas critical 
point [9] and, on the basis of our results, could be ex- 
tended also to the LLCP. However, this analysis is very 
expensive in terms of CPU time and goes beyond the goal 
of the present work. Nevertheless, the percolation analy- 
sis, presented in the next section, helps in understanding 
the physical reason for this large efficiency. 

The efficiency is a consequence of the fact that the av- 
erage size of Wolff's clusters changes with T and P in 



the same way as the average size of the regions of cor- 
related molecules (l9| . i.e. a Wolff's cluster statistically 
represents a region of correlated molecules. Moreover, 
the mean cluster size diverges at the critical point with 
the same exponent of the Potts magnetic susceptibility 
(lit . an d the clusters percolate at the critical point, as 
we will discuss in the next section. 



PERCOLATING CLUSTERS OF CORRELATED 
MOLECULES 

The efficiency of the Wolff's cluster algorithm is a con- 
sequence of the exact relation between the average size 
of the finite clusters and the average size of the regions 
of thermodynamically correlated molecules. The proof 
of this relation at any T derives straightforward from the 
proof for the case of Potts variables . This relation 
allows to identify the clusters built during the MC dy- 
namics with the correlated regions and emphasizes (i) the 
appearance of heterogeneities in the structural correla- 
tions 



22j , and (ii) the onset of percolation of the clusters 



of tetrahedrally ordered molecules at the LLCP [23|, as 
shown in Fig. [4] 



A systematic percolation analysis [LSj is beyond the 
goal of this report, however configurations such as those 
in Fig. 2] allow the following qualitative considerations. 
At T > Tc the average cluster size is much smaller than 
the system size. Hence, the structural correlations among 
the molecules extends only to short distances. This sug- 
gests that the correlation time of a local dynamics, such 
as Metropolis MC or molecular dynamics, would be short 
on average at this temperature and pressure. Neverthe- 
less, the system appears strongly heterogeneous with the 
coexistence of large and small clusters, suggesting that 
the distribution of correlation times evaluated among 
molecules at a given distance could be strongly heteroge- 
neous. The clusters appear mostly compact but with a 
fractal surface, suggesting that borders between clusters 
can rapidly change. 

At T ~ Tq there is one large cluster, in red on the right 
of the middle panel of Fig. [5J with a linear size compa- 
rable to the system linear extension and spanning in the 
vertical direction. The appearance of spanning clusters 
shows the onset of the percolation geometrical transition. 
At this state point the correlation time of local, such as 
Metropolis MC dynamics or molecular dynamics would 
be very slow as a consequence of the large extension of 
the structurally correlated region. On the other hand, 
the correlation time of the Wolff's cluster dynamics is 
short because it changes in one single MC step the state 
of all the molecules in clusters, some of them with very 
large size. Once the spanning cluster is formed, it breaks 
the symmetry of the system and a strong effective field 
acts on the molecules near its border to induce their reori- 
entation toward a tetrahedral configuration with respect 
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the molecules in the spanning cluster. 

As shown in Fig. 3, the spanning cluster appears as a 
fractal object, with holes of any size. The same large dis- 
tribution of sizes characterizes also the finite clusters in 
the system. The absence of a characteristic size for the 
clusters (or the holes of the spanning cluster) is the con- 
sequence of the fluctuations at any length-scale, typical 
of a critical point. 

At T < Tc the majority of the molecules belongs to 
a single percolating cluster that represents the network 
of tetrahedrally ordered molecules. All the other clus- 
ters are small, with a finite size that corresponds to the 
regions of correlated molecules. The presence of many 
small clusters gives a qualitative idea of the heterogene- 
ity of the dynamics at these temperatures. 



SUMMARY AND CONCLUSIONS 

We describe the numerical solution of mean field equa- 
tions and the implementation of the Wolff's cluster MC 
algorithm for a cell model for liquid water. The mean 
field approach allows us to estimate in an approximate 
way the phase diagram of the model at any state point 
predicting intriguing new results at very low T 21]. 

To explore the state points of interest for these predic- 
tions the use of standard simulations, such as molecular 
dynamics or Metropolis MC, is not effective due to the 
onset of the glassy dynamics [|| . To overcome this prob- 
lem and access the deeply supercooled region of liquid 
water, we adopt the Wolff's cluster MC algorithm. This 
method, indeed, allows to greatly accelerate the autocor- 
relation time of the system. Direct comparison of Wolff's 
dynamics with Metropolis dynamics in the vicinity of the 
liquid-liquid critical point shows a reduction of the auto- 
correlation time of a factor at least 10 4 . 

Furthermore, the analysis of the clusters generated 
during the Wolff's MC dynamics allows to emphasize how 
the regions of tetrahedrally ordered molecules build up on 
approaching the liquid-liquid critical point, giving rise to 
the backbone of the tetrahedral hydrogen bond network 
at the phase transition [23j . The coexistence of clusters 
of correlated molecules with sizes that change with the 
state point gives a rationale for the heterogeneous dy- 
namics observed in supercooled water [22| . 



FIG. 4: Three snapshots of the system with N = 100 x 100, 
showing the Wolff's clusters of correlated water molecules. 
For each molecule we show the states of the four arms and 
associate different colors to different arm's states. The state 
points are at pressure close to the critical value Pc (Pvo /e = 
0.72 ~ Pcvo/e) and T > T c (top panel, k B T/e = 0.0530), 
T ~ T c (middle panel, k B T/e = 0.0528), T < T c (bottom 
panel, ksT/e — 0.0520), showing the onset of the percolation 
at T ~ T c . 
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